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PACS 05 . 45 . -a - Nonlinear dynamics and Chaos 
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PACS 63 . 20 . Pw - Phonons in crystal lattices - Localized modes 

Abstract -We follow the dynamics of nonlinear waves in two-dimensional disordered lattices with 
tunable nonlinearity. In the absence of nonlinear terms Anderson localization traps the packet 
in space. For the nonlinear case a destruction of Anderson localization is found. The packet 
spreads subdiffusively, and its second moment grows in time asymptotically as t". We perform 
fine statistical averaging and test theoretical predictions for a. Along with a precise confirmation 
of the predictions in [Chemical Physics 375, 548 (2010)], we also find potentially long lasting 
intermediate deviations due to a growing number of surface resonances of the wave packet. 
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Introduction. — Anderson localization (AL) - the 
halt of wave propagation in random potentials due to ex- 
ponentially localized modes - was theoretically predicted 
over 50 years ago [I] and in the past decades, since ob- 
served within a variety of experiments, including optics 
[2}|5] and matter waves |6|7| . These two are of strong inter- 
est, in that AL can be strongly altered by nonlinear Kerr 
effects in disordered photonic lattices [8-10, or atomic 
Bose-Einstein condensate interactions in optical lattices 
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Research within nonlinear disordered media largely fo- 
cuses on wave packet evolution in one-dimensional (1-d) 
systems. Asymptotic subdiffusive spreading is observed. 
An extended debate of the characterizing power exponents 
appears to be clarified by the theoretical predic- 

for 1-d 
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tions and their numerical verifications in 
cases. Most studies focus on quartic nonlinearities which 
correspond to two-body interactions. Motivated greatly 

or at 



by experiment, e.g. in liquid crystal optics 28 29 
BEC-BCS crossovers in ultracold Fermi gases [30[|31| , one 
may also parametrize the nonlinearity exponent. This was 
done for the case of 1-d systems in 24 32 again with a 



confirmation of the theoretical prediction given in |25| . 
The innovation here is to extend to two-dimensional (2- 
d) disordered systems. For such lattices that are mul- 
tidimensional, disordered, and have variable nonlinearity 
exponents, spreading behaviors were broadly conjectured 
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within 25 . The aim of this letter is to (dis)proof the 



conjectures for asymptotic spreading in 2-d lattices with 
tunable nonlinearity. We will also investigate the case of 
small nonlinearity exponents at which the theory predicts 
an anomaly in the number of wave packet surface reso- 
nances which should grow with ongoing spreading. We 
will test the robustness of the theoretical predictions in 
this regime as well, where subdiffusion competes with fin- 
gering resonance instabilities. 

Generalized models. — The first model scrutinized 
is the generalized disordered nonlinear Schrodinger equa- 
tion (gDNLS), which in a discrete lattice reads 



nD = Y. 



\^rV + 



2/3|V'r 



E 



(^ArV'n+V'JV'n) 



(1) 

Here V'r are complex variables, where r = {x,y) denotes 
a 2-d square lattice vector of integer components with the 
nearest neighbor set of A/" C {r ± (1, 0), r ± (0, 1)}. The 
disorder appears in on-site energies Cr, which are uncor- 
related random values drawn uniformly from an interval 
[—W/2, W/2] parameterizing by the disorder strength W . 
The nonlinearity of strength j3 is generalized to a power 
CT > 0. This is best seen in the equations of motion, de- 
rived from ?/)r = dV.D/d{ii'*) as 
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The above set of dynamic equations conserves the total 
energy Hd, as well as the total norm S — J2r IV'rl • The 
1-d version of the gDNLS has been extensively studied: 
for cr = 2 it relates to recent experimental photonics [lO] 
and has been investigated numerically 23-27 . For a few 



integer values of a 1-d simulations were presented in 32 
Simulations with non-integer a were also performed 19 



on 



short time scales without focus on asymptotic spreading. 

The second model considered is the generalized 2-d 

Klein-Gordon (gKG) lattice, governed by the Hamiltonian 
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(3) 
where Uy. and p^ respectively are generalized coordinates 
and momenta on the lattice site r, with an energy den- 
sity of £r. The terms Cr are uncorrelated random val- 
ues drawn uniformly from an interval [1/2,3/2]. From 
Ur = —&HK/dur, the equations of motion read 

' - ' ■ ■ (4) 
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This set of dynamic equations conserve only the total en- 
ergy Hk = Sr^r- The 1-d version of the gKG has also 
been extensively studied, as it can be considered a model 
for dynamics of anharmonic optical lattice vibrations in 
molecular crystals [33]. The quartic 1-d case was heav- 



ily used in numerical investigations [23f 27 , and different 



values of a have also been addressed 34 



Several works 35 - 37 suggest an equivalence between 
the two models. Note, that in KG there is no nonlin- 
ear parameter /3. Rather the scalar value Hk acts as the 
nonlinearity control: by writing u^ as a plane wave in 
Eq.Q and applying slow modulation/rotating wave ap- 
proximations, Eq.M can be recovered under an approxi- 
mate condition. For quartic nonlinearity, this condition is 
f3S « 3WHk- It connects the KG initial parameters Hk 
and W to the total initial norm S and nonlinear param- 
eter /? of the corresponding quartic DNLS model. This 
condition can furthermore be generalized to any power a 
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(5) 
In this derivation, the absolute value in Eqs.(2J4| is ig- 
nored, since its inclusion was found only to yield minor 
higher order corrections to ao-. A similar result for the 
nonlinear shifts in energy was seen in [3 4|. 

Neglecting nonlinear terms both Eqs.(2J4| reduce to an 
eigenvalue problem, giving a set of exponentially local- 
ized eigenstates (denoted as normal modes, NM) with fre- 
quencies Ar in a spectrum of width Ad ^ 8 + W in the 
case of gDNLS. Linear reduction for the gKG is similar, 
but with squared frequencies w^ in a spectrum of width 
Ak — Ad/W = 1 -f 8/W. We will focus mainly on an- 
alytics of the gDNLS, since it is straightforward to adapt 
results for the gKG using Eqs.(l5|. 



Expected spreading regimes. — Consider the time- 
dependent normalized norm density distribution, z^ = 
\tjjr\ /S. The gKG counterpart is the normalized en- 
ergy density distribution, Zr = Sr/HK- Distributions 
are analyzed by means of the second moment, m2 = 
Sr I** ~ /^rl -^f where the density center is ^r = J2r ''•^r- 
The second moment quantifies the squared width of the 
packet, hence, its spreading. The participation number, 
P = 1/ J2r -^r I measures the number of effectively excited 
sites. Lastly, the packet sparseness is measured by the 
compactness index 23 , which for 2-d models is C = P/m2 ■ 



How then do these three measures behave for different pa- 
rameters? 

In the linear case, the participation number approxi- 
mates the spatial extension of a NM, with an average mea- 
sure over modes being the localization volume V. The de- 
pendence oiVonW is shown in the inset of FigjT] where 
squares [diamonds] are for the linear version of Eq.(IT]) 
[Eq.(l3])], the gray cloud is the overall standard deviation, 
and the solid line is a best fit of y ~ W through all 
points. Similar curves also appear in 38 39 . The average 



frequency spacing of NMs within the localization volume 
is d = Ad /v. The two linear frequency scales d and Ad 
are expected to contribute to the details of packet spread- 
ing. Nonlinearity also introduces an additional frequency 
scale - the nonlinear shift of a single oscillator, propor- 
tional to /3/9°'/^ for the gDNLS, where p is the average 
norm density of a packet. From these three unique fre- 
quency scales, dynamical regimes of packet spreading were 



presented in 1-d quartic systems 23-27 for a variety of 



different initial parameters. Both the initial norm/energy 
density of a packet and its typical size were suggested 
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as the major control parameters for the dynamics 
at given W and a. Under strong enough nonlinearity, 
a fraction of a wave packet (or even the whole packet) 
exhibits self-trapping 26 27 40]. For weaker nonlinearity 



(such that self-trapping is avoided), two possible dynami- 
cal outcomes are predicted. The packet spreads in an in- 
termediate regime of strong chaos with subsequent dynam- 
ical crossover into an asymptotic regime of weak chaos, or 
spreading starts directly in the weak chaos regime. 

A straightforward generalization of these expected 
regimes of packet spreading, with an initial packet norm 
density p and size L < V, was proposed [25^ as 



pp^^^L/vy^^ < d 

Pp''/'^{L/Vf/^ > d 
/3p"/2 > Ad 



weak chaos, 
strong chaos, 
self-trapping. 



(6) 



The spreading mechanism is thought to be an incoherent 
energy transfer between NMs inside the packet to nearby 
exterior NMs 
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The derivation of possible dy- 
namical outcomes is based on the resonance probability V 
of this transfer; this in turn depends on the norm/energy 
density of the packet as Vi/Sp"/'^) « 1-exp {-/Sp'^/'^d-^). 
In the regime of weak chaos, only a small fraction of in- 
terior modes resonantly interact and V « f3p''/^d~^. In 



p-2 



Universal subdiffusion of nonlinear waves in two dimensions with disorder 



the regime of strong chaos, nearly all modes interact, 
i.e. V ^ 1. The diffusion rate D is conjectured to be 
D ~ /32p'^(-p(/3p'^/2))2 wiiich together with m2 -- p"\ 
leads to power laws for spreading in 2-d 



m2,P^f' 



a 




weak chaos, 
strong chaos. 



(7) 



Note, since the packet norm/energy density decreases in 
time and eventually the condition for strong chaos will be 
no longer satisfied - spreading will cross into the regime 
of weak chaos. Nevertheless, the duration of strong chaos 
regime can be greatly prolonged (over multiple orders of 
magnitude), so much that the crossover occurs at infeasi- 
ble computation times. 

In the regime of weak chaos, the number of resonances 
in the packet volume Njiy and on its surface Nrs are 



estimated 25 as 
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According to the above equation, for the 2-d case there 
are critical values of nonlinearity power: the number of 
volume resonances will grow for any ct < 2, likewise cr < 1 
for surface resonances. We therefore expect these critical 
values may manifest unusual effects in the course of packet 
spreading (we shall return to this point in the last section). 

Numerical Simulations. — Our following numerics 
present only gKG results, for two reasons. Firstly, the 
presence of a corrector scheme for gKG (e.g. Appendix of 
[24| ) allows two magnitude orders greater in integration, 
at the same conservations and integration speeds. Sec- 
ondly, the gKG requires only a single conservation. Lastly, 
all prior simulations of both models [23}j27l show similar 
qualitative results in a wide range of energy and disorder. 

To test the analytical prediction of Eq.(l7]), we first set 
the localization volume to y ~ 34 for W = 10. This will 
stay the same for all presented numerics. For an initial 
single-site excitation {L = 1) of energy £, regime bound- 
aries from Eq.Q can be easily mapped into a £{(t) form 
using Eq.Q. This effectively gives a parameter space, 
shown in the main of Fig [I] The two dashed lines show 
the boundary deviation, due to variations in V (see inset). 

We then select energy density £ so that our initial state 
is set as pr — v2£ for r = (-^/2, N/2) and zero else- 
where, Mr = everywhere. This state is numerically inte- 
grated by Eq.Q using SABA-class symplectic integrators 
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under a time-step of 10~^ to maximums of 10^~®, all 
the while maintaining a relative energy conservation up to 
10~^. For each set of parameters, we calculate the three 
measures m2,P,C- The typical lattice size is 200 x 200 
sites. We then average over 400 disorder realizations. For 
each parameter set, we determine the spreading power-law 
exponent as local derivative a(t) = d (log j^q 7712) /dlogj^p^ 
(see, e.g. [26]). 



Fig. 1: Main: parameter space of the nonlinearity power a and 
the energy £ of a single-site excitation (L = 1). Dashed lines 
show the variation of the boundary obtained from variation 
of V (see inset). For our numerics, the different symbols 
correspond to points (a, £) used. The various behaviors dis- 
cussed in further sections of the text: "*" for (2, 0.3), (2, 2.0), 
"-h" for (0.5, 0.00001), (0.7, 0.0005), (1.0, 0.006), "o" for 
(1.3, 0.025), (1.5, 0.04), and "x" for (0.5, 0.005), (0.7,0.03). 
Inset: the dependence of localization volume V on W. Squares 
[diamonds] are for the linear version of Eq. (II]) [Eq. S]. The 
gray region denotes an overall standard deviation. The solid 
line is a best fit to \^ ~ W^ . 



Subdiffusion in 2-d lattices with a — 2: weak chaos 
and self-trapping. In this subsection, the numerics for 
quadratic nonlinearity ("*" in Figjl]) are discussed. With 
single-site excitations, the strong chaos regime is unreach- 
able for CT > 2. In FigJ2] the weak chaos regime {£ = 0.3, 
lower "*" in FigfTl) displays spreading, with TO2 and P 
growing (blue curves) to reach an asymptotic spreading 
following the predicted power law of Eq.Q, seen by the 
saturation a ~ 0.2055 and verifying our theoretical pre- 
diction of 1/5 for weak chaos. This ought be compared 
to an earlier work [43| , in which the asymptotic law was 
hypothesized as 1/4, but only via integration to 10^. In- 
deed, at this time, we also approximately obtain the same 
value (a ~ 0.234); however upon integrating further, our 
expectation of a = 1/5 is revealed by the saturation. In 
this regime, the asymptotic compactness index is C — 2.36 
(seen in the blue curve of FigJ2]) , meaning that the packet 
spreads, yet remains largely thermalized (C « 3). A trend 
of C — >■ indicates either very sparse packets or partial self- 
trapping [26]. Large values of £ satisfy the self-trapping 
regime. This behavior can be seen in the red curves of 
Figl2| {£ — 2.0, upper "*" in Figjl]). A large portion of 
energy remains trapped on the initially excited site, while 
a much smaller portion subdiffuses. We thusly observe an 
increase of 7712 (smaller portion subdiffusing) , while P re- 
mains largely unaffected (large self-trapped portion). The 
compactness index approaches zero - a very good indica- 
tion of self-trapping. Lastly, in the inset of Fig|2] we show 
normalized density distributions (averaged first over real- 
ization, then over the y-coordinate) for both regimes. The 
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Fig. 2: (Color online.) Numerics for "*" in Figll] The 
parameters {<y,£) ~ (2,0.3), (2, 2.0) correspond to the weak 
chaos ((b)lue) and self-trapping ((r)ed). Left column: sec- 
ond moment (upper) and its power-law exponent a (lower). 
The dashed line is our theoretical expectation for weak chaos 
a = 1/(1 + 2(t) = 0.20. Right column: participation number 
(upper) and compactness index (lower). In both columns of 
the upper row the lighter clouds correspond to a standard de- 
viation. Inset: Normalized density distributions at f = 10*, 
averaged first over realization and then over the y-coordinate. 
Self-trapping is clearly seen. 



self-trapping regime reveals a characteristically trapped 
portion in the center. 

Subdiffusion in 2-d lattices with 1 < a < 2: weak chaos. 

Numerical findings for "o" in Figjl] with parameters 
{a,£) = (1.3, 0.025), (1.5, 0.04) are shown in Figjs) Par- 
ticularly in the lower left panel, the exponent a is shown 
to also asymptotically saturate. The two "I-bars" do not 
show error per se, rather their lower/upper bounds dic- 
tate the weak/strong chaos expectations for a from Eq.(l7]). 
The asymptotic saturations are approximately a ~ 0.300 
for a — 1.3 and a ~ 0.257 for a — 1.5, which is quite close 
to their respective weak chaos expectations of 0.278 and 
0.250. Additionally, compactness values of C — 2.44,2.10 
at i = lO'' remain fairly thermalized. Therefore, these 
two points fully follow our expected theory of weak chaos 
spreading. 

Subdiffusion in 2-d lattices with a < 1: strong chaos. 

Moving to the left in FigjT] we cross the theoreti- 
cal division between strong and weak chaos (two rep- 
resentative points are shown by "x"). Again, we per- 
form our numerics for these two points. Note, that red 
curves are for {a,£) = (0.5,0.005) and blue curves are 
for {a,£) = (0.7,0.03). Similarly, the I-bar bounds give 
weak chaos (upper) and strong chaos (lower) expectations 
for the respective values of a. The asymptotic satura- 
tions for these two representative points are approximately 
a ~ 0.669 for a = 0.5 and a ~ 0.571 for a = 0.7, which is 
quite close to their respective strong chaos expectations of 
0.667 and 0.589. Additionally, the compactness index flue- 



Fig. 3: (Color online.) Numerics for "o" in FigfT] The parame- 
ters {(J,£) = (1.3,0.025), (1.5,0.04) are colored respectively as 
(r)ed and (b)lue. Left column: the second moment (upper) 
and its power-law exponent a (lower). The I-bar bounds de- 
note the theoretical expectations from Eq.([7| for weak chaos 
(lower bound) and strong chaos (upper bound). Right column: 
participation number (upper) and compactness index (lower). 
In both columns of the upper row the lighter clouds correspond 
to a standard deviation. 



tuates due to a slight asymptotic slope change in the par- 
ticipation numbers, nevertheless, it remains nearly ther- 
malized about 1.7 for t — 10^. Therefore, these two points 
follow our expected theory of strong chaos spreading. 




Fig. 4: (Color online.) Numerics for "x" in FigfT] The param- 
eters {(J,£) = (0.5,0.005), (0.7, 0.03) are colored respectively 
as (r)ed and (b)lue. Left column: the second moment (up- 
per) and its power-law exponent a (lower). The I-bars denote 
the theoretical expectations from Eq.([7| for weak chaos (lower 
bound) and strong chaos (upper bound). Right column: par- 
ticipation number (upper) and compactness index (lower). In 
both columns of the upper row, the lighter clouds correspond 
to a standard deviation. 

Subdiffusion in 2-d lattices with a < 1: intermediate 
behaviors. Moving down in FigfT] we cross back to the 
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theoretical weak chaos regime (three representative points 
are given by "+" ) . Performed numerics are shown in FigJ5| 
where red curves are for (tr, £) = (0.5,0.00001), green 
curves are for {a, £) = (0.7, 0.0005), and blue curves are for 
(cr, f) — (1.0,0.006). These points exhibit behavior novel 
from that was seen in the previous sections. Namely, in the 
lower left panel the exponent power a reaches asymptotic 
values of 0.586, 0.494, 0.375 respectively for red, green, and 
blue curves. That is a clear tendency toward saturations 
resting in neither regime. Rather, these values lay be- 




Fig. 5: (Color online.) Numerics for "+" in Figfl] The pa- 
rameters {a,£) = (0.5, 0.00001), (0.7, 0.0005), (1.0, 0.006) are 
colored respectively as (r)ed, (g)reen, and (b)lue. Left col- 
umn: the second moment (upper) and its power-law expo- 
nent Q (lower). Similarly, I- bar bounds denote the theoreti- 
cal expectations from Eq.|[7| for weak chaos (lower bound) and 
strong chaos (upper bound) . Right column: participation num- 
ber(upper) and compactness index (lower). In both columns 
of the upper row, the lighter clouds correspond to a standard 
deviation. 

tween the two bounds from the expected regimes. Similar 
behavior was also hinted about cr = 1.0 for 1-d, particu- 
Returning to our argument of resonance 



larly Fig.5 of 34 



probability, this suggests both the weak/strong limits are 
invalid - the value of V must be explicitly found. More 
than just a few modes contribute, but certainly not enough 
to yield strong chaos regime. 

Dimensional Analysis. According to Eq.([8|, the con- 
cept of resonance probability may be viewed in the light 
of competition between surface growth versus volume 
growth. Surface resonances more easily lead to density 
leakage into modes exterior to the packet. This process 
in turn increases the packet's perimeter, therefore yield- 
ing more surface resonances. Packets may thusly develop 
finger structures or fragment, perhaps leading to a fractal- 
like structure. 

As a first response into this, we consider the normalized 
densities z^ — St/Hk at i = 10®, where an asymptotic 
regime is reached (cf. FigsJ2|5|. In Fig|6J the largest con- 
tour path for Zr > 10"^ is shown, corresponding roughly 



to the expanding packet's surface. The left panel is for 
weak chaos with ct > 1, the middle panel is for weak 
chaos with <t < 1, and the right panel is for strong chaos. 
Overlaid in black is the localized linear packet surface. 
A comparison of the coefficient a in Figsj3][5] can be ob- 
served. Wave packet in the left panel (a corresponds to 
weak chaos asymptotic) spreads less than presented in the 
middle panel (a in between the two limits) , which spreads 
less than the right panel {a corresponds to strong chaos 
limit). However, the boundary shape itself provides no 
further evidence: no one boundary appears to be more 
fragmented or fingered than the others. 



Regime 


{a,£) 


(Df) 


Linear 


N/A 


1.498 ± 0.045 


Weak Chaos, a > 1 


(2.0,0.3) 
(1.5,0.04) 
(1.3,0.025) 


1.621 ±0.019 
1.601 ±0.021 
1.623 ±0.017 


Weak Chaos, ct < 1 


(1.0,0.006) 
(0.7,0.0005) 
(0.5,0.00001) 


1.667 ±0.015 
1.723 ±0.011 
1.740 ±0.007 


Strong Chaos 


(0.7,0.030) 
(0.5,0.005) 


1.734 ±0.009 
1.730 ±0.007 



Table 1: Box-Counting Dimension for the different regimes, 
Zr > 10-^ 



Therefore, we turn to a dimensional analysis of the 
contours. Binarizing z^ at the the threshold > 10"^, a 



box-counting algorithm 44 is performed to extract the 



Minkowski-Bouligand dimension Df of the surface. The 
results were averaged over 100 different realizations and 
presented in Table [l] This further suggests there are no 
clear and distinct fragmentation/fingering structures that 
might separate (in a geometric sense) the weak chaos at 
CT < 1 from the other two regimes. 



Conclusion. — We have investigated the spreading of 
a single-site excitation under a variable power nonlinearity 
within a 2-d disordered lattice, in particular for the Klein- 
Gordon case of Eq.([3]). For such a system, we numerically 
confirm the second moment behavior of Eq.Q, as first 
hypothesized in |25| . In particular, we verify existence of 
both a weak chaos and a strong chaos regime. In addition, 
an intermediate regime for ct < 1 is observed, with spread- 
ing behavior between the two limits of strong and weak 
chaos. We have performed an analysis of the wavepacket 
geometries, but so far, strong fragmentation/fingering ev- 
idence in this regime, compared to the other regimes, is 
eluding. Possible future avenues along these lines may in- 
clude lucunarity and density-density correlation measures. 
The behavior between the two regimes certainly remains 
open for future exploration, as well as pushing numerically 
the DNLS to achieve similar observations. 
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Fig. 6: Largest contour for Zr > 10~^ at t = 10®. The black area in the center is the contour for the linear case, while the panels 
respectively correspond to (from left to right): weak chaos with {(JjS) — (1.5,0.04), weak chaos with {cr,£) — (0.5,0.00001), 
and strong chaos {<J,£) ~ (0.5,0.005). 
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